# Time Series

# 生成1-30的整数，频率是12也就是月度数据，从2011年3月开始
a <- ts(1:30, frequency = 12, start = c(2011, 3))
print(a)

str(a)
attributes(a)

# AirPassengers

plot(AirPassengers)

apts <- ts(AirPassengers, frequency = 12)

# 使用函数decompose()分解时间序列
f <- decompose(apts)

# 季度数据
f$figure
plot(f$figure, type = "b", xaxt = "n", xlab = "")

# 使用当前的时间区域给月份赋予名称
monthNames <- months(ISOdate(2011, 1:12, 1), abbreviate = F)

# 使用月份名称标记x轴
# side=1表示设置x轴，at指的是范围从10-12，las表示分割的单位刻度为2
axis(1, at = 1:12, labels = monthNames, las = 2)
plot(f)

# Predict
# 参数order由（p,d,q）组成，p=1指的是自回归项数为1，list内容是季节seasonal参数
fit <- arima(AirPassengers, order = c(1, 0, 0), 
             list(order = c(2, 1, 0), period = 12))
# 预测未来24个月的数据
fore <- predict(fit, n.ahead = 24)
# 95%的置信水平下的误差范围（U,L）
u95 <- fore$pred + 2 * fore$se
l95 <- fore$pred - 2 * fore$se
# col=c(1,2,4,4)表示线的颜色分别为黑色，红色，蓝色，蓝色
# lty=c(1,1,2,2)中的1，2指连接点的先分别为实线和虚线
ts.plot(AirPassengers, fore$pred, u95, l95, 
        col = c(1, 2, 4, 4), lty = c(1, 1, 2, 2))
legend("topleft", c("Actual", "Forecast", "Error Bounds (95% CI)"), 
       col = c(1, 2, 4), lty = c(1, 1, 2))

# Dynamic Time programming

library(dtw)
idx <- seq(0, 2 * pi, len = 100)


a <- sin(idx) + runif(100) / 10
b <- cos(idx)

align <- dtw(a, b, step = asymmetricP1, keep = T)
dtwPlotTwoWay(align)

# 

sc <- read.table("synthetic_control.data", header = F, sep = "")

# 显示每一类数据的第一个样本观测值
idx <- c(1, 101, 201, 301, 401, 501)
sample1 <- t(sc[idx, ])

set.seed(6218)
n <- 10
s <- sample(1:100, n)
idx <- c(s, 100 + s, 200 + s, 300 + s, 400 + s, 500 + s)
sample2 <- sc[idx, ]
observedLabels <- rep(1:6, each = n)

# 使用欧式距离层次聚类
hc <- hclust(dist(sample2), method = "average")
plot(hc, labels = observedLabels, main = "")

# 将聚类树划分为6类
rect.hclust(hc, k = 6)
memb <- cutree(hc, k = 6)
table(observedLabels, memb)

## DTW distance

library(dtw)

distMatrix <- dist(sample2, method = "DTW")
hc <- hclust(distMatrix, method = "average")
plot(hc, labels = observedLabels, main = "")
rect.hclust(hc, k = 6)
memb <- cutree(hc, k = 6)
table(observedLabels, memb)


# 给原始的数据集加入分类标签classId
classId <- rep(as.character(1:6), each = 100)
newSc <- data.frame(cbind(classId, sc))

library(party)
ct <- ctree(classId ~ ., data = newSc, 
            controls = ctree_control(minsplit = 30, minbucket = 10, maxdepth = 5))
pClassId <- predict(ct)
table(classId, pClassId)

# 计算分类的准确率
(sum(classId == pClassId)) / nrow(sc)
plot(ct, ip_args = list(pval = F), ep_args = list(digits = 0))


# DWT

library(party)
library(wavelets)

wtData <- NULL

# 遍历所有时间序列
for (i in 1:nrow(sc)) {
  a <- t(sc[i, ])
  wt <- dwt(X = a, filter = "haar", boundary = "periodic")
  wtData <- rbind(wtData, unlist(c(wt@W, wt@V[[wt@level]])))
}

wtData <- as.data.frame(wtData)
wtSc <- data.frame(cbind(classId, wtData))

# 使用DWT建立一个决策树，control参数是对决策树形状大小的限制
ct <- ctree(classId ~ ., data = wtSc, 
            controls = ctree_control(minsplit = 30, minbucket = 10, maxdepth = 5))
pClassId <- predict(ct)

table(classId, pClassId)


# KNN

k <- 20
# 通过在第501个时间序列中添加噪声来创建一个新的数据集
newTS <- sc[501, ] + runif(100) * 15
# 使用‘DTW’方法计算新数据集与原始数据集之间的距离
distances <- dist(newTS, sc, method = "DTW")
# 给距离升序排列
s <- sort(as.vector(distances), index.return = TRUE)
# s$ix[1:k]是排行在前20的距离，表哥输出k个最近邻居所属的类
table(classId[s$ix[1:k]])
